This example explains how to use the OpenPNM.Utilies.IO.Statoil class to import a network produced by the Maximal Ball network extraction code developed by Martin Blunt's group at Imperial College London. The code is available from him upon request, but they offer a small library of pre-extracted networks on their website.
In [1]:
import warnings
import scipy as sp
import numpy as np
import openpnm as op
np.set_printoptions(precision=4)
np.random.seed(10)
%matplotlib inline
The following assumes that the folder containing the 'dat' files is in a directory called 'fixtures' in the same directory as this script. You can also enter a full path to the files.
In [2]:
from pathlib import Path
path = Path('../fixtures/ICL-Sandstone(Berea)/')
project = op.io.Statoil.load(path=path, prefix='Berea')
pn = project.network
pn.name = 'berea'
This import class extracts all the information contained in the 'Statoil' files, such as sizes, locations and connectivity. Note that the io classes return a project object, and the network itself can be accessed using the network attribute. The following printout display which information was contained in the file:
In [3]:
print(pn)
At this point, the network can be visualized in Paraview. A suitable '.vtp' file can be created with:
In [4]:
project.export_data(filename='imported_statoil')
The resulting network is shown below:
Although it's not clear in the network image, there are a number of isolated and disconnected pores in the network. These are either naturally part of the sandstone, or artifacts of the Maximal Ball algorithm. In any event, these must be removed before proceeding since they cause problems for the matrix solver. The easiest way to find these is to use the check_network_health function on the network object. This will return a dictionary with several key attributes, including a list of which pores are isolated. These can then be trimmed using the trim function in the topotools module.
In [5]:
print('Number of pores before trimming: ', pn.Np)
h = pn.check_network_health()
op.topotools.trim(network=pn, pores=h['trim_pores'])
print('Number of pores after trimming: ', pn.Np)
When importing Statoil networks, OpenPNM must perform some 'optimizations' to make the network compatible. The main problem is that the original network contains a large number of throats connecting actual internal pores to fictitious 'reservoir' pores. OpenPNM strips away all these throats since 'headless throats' break the graph theory representation. OpenPNM then labels the real internal pores as either 'inlet' or 'outlet' if they were connected to one of these fictitious reservoirs.
It is fairly simple to add a new pores to each end of the domain and stitch tehm to the internal pores labelled 'inlet' and 'outlet', but this introduces a subsequent complication that the new pores don't have any geometry properties. For this example, we will not add boundary pores, but just the pores on the inlet and outlet faces.
In OpenPNM 2+ it is optional to define Geometry and Physics objects (These are really only necessary for simulations with diverse geometrical properties in different regions, resulting in different physical processes in each region, such as multiscale networks for instance). It is still necessary to define Phase objects:
In [6]:
water = op.phases.Water(network=pn)
In [7]:
water.add_model(propname='throat.hydraulic_conductance',
model=op.models.physics.hydraulic_conductance.valvatne_blunt)
Recall that boundary pores and throats had no geometrical properties associated with them, so the hydraulic conductances of boundary throats will be undefined (filled with NaNs):
In [8]:
print(water['throat.hydraulic_conductance'])
In [9]:
flow = op.algorithms.StokesFlow(network=pn, phase=water)
flow.set_value_BC(pores=pn.pores('inlets'), values=200000)
flow.set_value_BC(pores=pn.pores('outlets'), values=100000)
flow.run()
The resulting pressure field can be visualized in Paraview, giving the following:
In [10]:
# Get the average value of the fluid viscosity
mu = np.mean(water['pore.viscosity'])
# Specify a pressure difference (in Pa)
delta_P = 100000
# Using the rate method of the StokesFlow algorithm
Q = np.absolute(flow.rate(pores=pn.pores('inlets')))
# Because we know the inlets and outlets are at x=0 and x=X
Lx = np.amax(pn['pore.coords'][:, 0]) - np.amin(pn['pore.coords'][:, 0])
A = Lx*Lx # Since the network is cubic Lx = Ly = Lz
K = Q*mu*Lx/(delta_P*A)
print(K)